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Abstract —Targets that generate multiple measurements at a 
given instant in time are commonly known as extended targets. 
These present a challenge for many tracking algorithms, as they 
violate one of the key assumptions of the standard measurement 
model. In this paper, a new algorithm is proposed for tracking 
multiple extended targets in clutter, that is capable of estimating 
the number of targets, as well the trajectories of their states, 
comprising the kinematics, measurement rates and extents. The 
proposed technique is based on modelling the multi-target state 
as a generalised labelled multi-Bernoulli (GLMB) random finite 
set (RFS), within which the extended targets are modelled 
using gamma Gaussian inverse Wishart (GGIW) distributions. A 
cheaper variant of the algorithm is also proposed, based on the 
labelled multi-Bernoulli (LMB) filter. The proposed GLMB/LMB- 
based algorithms are compared with an extended target version 
of the cardinaUsed probability hypothesis density (CPHD) filter, 
and simulation results show that the (G)LMB has improved 
estimation and tracking performance. 

Index Terms —Random finite sets, finite set statistics, GLMB 
filter, extended targets, multi-target tracking. Inverse Wishart, 
CPHD filter 


I. Introduction 

Multi-target tracking is the process of estimating the num¬ 
ber of targets and their states, based upon imperfect sensor 
measurements that are typically corrupted by noise, missed 
detections, and false alarms. The main challenge is to filter out 
these three effects in order to gain accurate estimates of the 
true target states. A Bayesian approach to this type of problem 
requires models to describe how the measurements are related 
to the underlying target states. Most traditional trackers use the 
so-called standard measurement model. This is also known as 
a point target model, since it is based on the assumption that 
each target produces at most one measurement at a given time, 
and that each measurement originates from at most one target. 
This model simplifies the development of multi-target trackers, 
but in practice it is often an unrealistic representation of the 
true measurement process. 

More realistic measurement processes can be handled using 
more sophisticated non-standard measurement models, which 
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may relax the aforementioned assumptions, usually at the ex¬ 
pense of increased computation. One example of this is when 
a group of targets produces a single measurement, known as 
an unresolved target (or merged measurement) model |T]. This 
model is useful when dealing with low-resolution sensors that 
cannot generate seperate detections for closely spaced targets. 
On the other hand, higher resolution sensors may produce 
multiple measurements per target on any given scan. Such 
cases require the use of an extended target model Q, which 
is the subject of this paper. 

Extended target measurement models typically require two 
components; a model for the number of measurements gener¬ 
ated by each target, and a model for their spatial distribution. 
These depend strongly on both the sensor characteristics and 
the type of targets being tracked. For example, in radar 
tracking, some targets may generate many separate detections, 
by virtue of the fact that they possess many scatter points. 
However, other targets may reflect most of the energy away 
from the receiver, leading to very few detections, or none at 
all. In general, when a target is far enough away from the 
sensor, its detections can often be characterised as a cluster of 
points exhibiting no discernable geometric structure. In such 
cases, the number of measurements is usually modelled using 
a Poisson distribution, see for example Q and 0. 

Even in the absence of a specific target structure, it is still 
possible to estimate the size and shape of a target, known as the 
target extent. This can be achieved by assuming some general 
parameteric shape for the extent, for which the parameters are 
estimated based on the spatial arrangement of the observations. 
An approach that assumes an elliptical extent was proposed in 
E’ which used a multivariate Gaussian, parameterised by a 
random covariance matrix with an inverse Wishart distribution. 
This was termed a Gaussian inverse Wishart (GIW), and this 
method enables the target extent to be estimated on-line, in¬ 
stead of requiring prior specification. Further applications and 
improvements have appeared in 0-0- Alternative methods 
for estimating target extent have also been proposed, see for 
example |8|- ID- 

The GIW method has been applied using multi-target filters 
based on the random finite set (RFS) framework. A proba¬ 
bility hypothesis density (PHD) filter, which was originally 
developed by Mahler in GD for the point target model, 
was proposed for extended multi-target filtering in |12|. An 
implementation of this filter based on the GIW model (GIW- 
PHD filter) was developed in GD- The cardinalised PHD 
(CPHD) filter | |27) is a generalisation of the PHD filter, which 
models the multi-target state as an i.i.d cluster RFS instead of 
a Poisson RFS. This was applied to extended targets in G3’ 
which also incorporated a modification to the GIW approach 
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HD. enabling the estimation of target measurement rates. This 
method treats the rate parameter of the Poisson pdf (which 
characterises the number of measurements generated by a 
target) as a random variable, whose distribution is modelled as 
a gamma pdf. This algorithm was called the gamma Gaussian 
inverse Wishart CPHD (GGIW-CPHD) filter. Extended target 
PHD and CPHD filters have also been presented in p2) , p^ . 

The advantage of the (C)PHD filters is that they reduce the 
computational cost of the Bayes multi-target filter, but to do 
so, some significant approximations are made. While these 
approximations avoid explicit data association, it means that 
the filters do not produce target tracks, and the PHD filter 
can produce highly uncertain estimates of the target number 
due to the Poisson cardinality assumption |[TT),|T8). Another 
limitation affecting the CPHD filter is the so-called ‘spooky’ 
effect Gl), which means a target misdetection may cause a 
false estimate to spontaneously appear in a different part of 
the state space. A Bernoulli filter for extended targets was 
proposed in pO) , which does not suffer from these issues, 
however, it is limited to at most a single target in clutter. 

A recently proposed algorithm that addresses these lim¬ 
itations is called the generalised labelled multi-Bernoulli 
(GLMB) filter pT| , | |22) . This algorithm has been shown to 
outperform both the PHD and CPHD filters, with the added 
advantage of producing labelled track estimates, albeit with 
a higher computational cost. An approximate but computa¬ 
tionally cheaper version of this filter was proposed in p3) , 
called the labelled multi-Bernoulli filter (LMB). Also, the 
first GLMB filter for a non-standard measurement model was 
developed in ||Tj, which used a model that includes merged 
measurements, a problem which can be viewed as the dual of 
the extended target tracking problem. 

In this paper, we develop a GLMB filter for extended multi¬ 
target tracking based on the GGIW model. The resulting algo¬ 
rithm (GGIW-GLMB) is capable of estimating the kinematics 
and extents of multiple extended targets in clutter, with the 
advantage of producing full target tracks. Preliminary results 
on this work have been presented in p^ , which we build 
upon in the following ways. Lirstly, we provide a complete 
derivation of the extended target likelihood function used 
by the filter. Second, we have implemented and tested a 
computationally cheaper version of the algorithm, called the 
GGIW labelled multi-Bernoulli (GGIW-LMB) filter. Third, we 
have improved the utility of the filter by incorporating an 
adaptive target birth model, allowing new targets to appear 
from anywhere in the state space. Linally, we have applied 
the algorithms to a real-world data set obtained from a lidar 
sensor used in autonomous vehicle applications. 

The paper is organised as follows. Section contains a 
brief background on the GLMB, and in Section [III] we adapt 
it to extended multi-target tracking, by proposing an extended 
target likelihood and deriving the associated GGIW-GLMB 


and GGIW-LMB filters. Section IV presents details relating 


to the implementation of these algorithms. Section [V| contains 
simulation results comparing the performance of the GGIW- 
(G)LMB with the GGIW-CPHD filter, and in Section we 
demonstrate an application to real-world measurement data. 


H. Background; Tracking with Labelled Random 
Linite Sets 

The essence of the RLS approach to multi-target tracking is 
the modelling of the multi-target states and measurements as 
finite set-valued random variables, or RLSs. Until recently, 
estimation algorithms derived using this framework have 
been based on the use of unlabelled random finite sets, as 
demonstrated by the PHD fTT] , CPHD p7) and multi-target 
multi-Bernoulli (MeMBer) |28| filters. A key reason for the 
popularity of these approaches is that they do not require 
explicit data association. However, their main disadvantage is 
that they only provide sets of unlabelled point estimates at each 
time step, so in applications that require target trajectories, 
tracks must be formed via additional post-processing. To 
address this problem, the concept of labelled random finite 
sets was proposed lH), which involves assigning a distinct 
label to each target, such that each target’s trajectory can be 
identified without the need for post-processing. 

In pT| , a class of labelled RLS called generalised labelled 
multi-Bernoulli (GLMB) was proposed, and based on this 
formulation, an algorithm for solving the multi-target tracking 
problem under the standard point-target likelihood model was 
developed. In the remainder of this section we briefly review 
some of the key points of this technique, and in Section III 
we adapt this method to tracking multiple extended targets. 

We begin by introducing some notation and definitions relat¬ 
ing to labelled random finite sets. The multi-object exponential 
of a real valued function h raised to a set X is defined as 
where ft,® = 1, and the elements 
of X may be of any type such as scalars, vectors, or sets, 
provided that the function ft(-) takes an argument of that type. 
The generalised Kronecker delta function, and the set inclusion 
function are respectively defined as 


(X) = 


if X = Y 
otherwise 


If (X) = 


1, if X c r 

0, otherwise 


where again, X and Y may be of any type, such as scalars, 
vectors, or sets. In general, we adopt the notational convention 
that labelled sets are expressed in bold upper case (X), 
unlabelled sets in regular upper case (X), labelled vectors 
in bold lower case (x), and unlabelled vectors or scalars in 
regular lower case (x). 

Definition 1. A labelled RFS X with state space X and 
discrete label space L, is an RLS on X x L, such that the 
labels within each realisation are always distinct. That is, if 
£(X) is the set of unique labels in X, and we define the 
distinct label indicator function as 


A(X) = 


1, if|£(X)| = |X| 


^0, if|£(X)|^|X| 
then a labelled RLS X always satisfies A (X) = 1. 


( 1 ) 


Definition 2. A labelled multi-Bernoulli (LMB) RLS is a 
labelled RLS with state space X and discrete label space L, 
which is distributed according to 


Linally, we make some concluding remarks in Section VII 


7r(X) = A(X)u;(/:(X))[p(.)f , 


( 2 ) 
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where 

»(L)=n(i-^<‘>)n^ 

ieh l&L 

pix,l) (x ), 

in which and p^^'> (•) are the existence probability and 
probability density corresponding to label I G L. An LMB 
distribution is abbreviated using the notation tt (JC) — 

{(r(0,p(0)}^^^. 

Definition 3. A generalised labelled multi-Bernoulli (GLMB) 
RFS is a labelled RFS with state space X and discrete label 
space L, which is distributed according to 


(3) 

(4) 


TrfX) 


cGC 




where C is a discrete index set, and (L) and p^‘^'1 {x,l) 
satisfy 


^ ^ (L) = 1, / p^^'^ {x,l)dx = l. (6) 

LCLcGC 


In Bayesian multi-target tracking, the goal at time k is to 
estimate a finite set of labelled states C X x L, called 
the multi-target state, based on finite sets of incoming multi¬ 
target observations Zk C Z. We model as a labelled 
random finite set, and Zk as an unlabelled random finite set. 
A principled mathematical framework for working with RFSs 
is called finite set statistics (FISST) | |25l , the cornerstone of 
which is a notion of multi-target density/integration that is 
consistent with point process theory p6) . 

The multi-target state at each time k is distributed according 
to a multi-target density tt^ (- 1 X 1 :^), where Zi,k is an array 
of finite sets of measurements received up to time k. The 
multi-target density is recursively propagated in time via a 
multi-target prediction and update as follows. 

The multi-target prediction to time k is given by the 
Chapman-Kolmogorov equation 


tfk\k-i {Xk\Zi.k-i) (7) 

= I fk\k-i{Xk\X)7Tk-i{X\Zi.,k-i)5X, 

where fk\k-i {Xk\X) is the multi-target transition kernel 
from time fc — 1 to time k, and the integral is the set integral, 
defined by for any function / that takes (X x L), the 
collection of all finite subsets of X x L, to the real line. 


/ (X) sx = 


oo 

E 

2 = 0 




( 8 ) 

At time k, a set of observations Zk is received, which is 
modelled by a multi-target likelihood function gk{Zk\Xk)- 
The multi-target posterior at time k is given by Bayes rule 


TTk (-^fc|X^l:fc) — 


gk {Zk\Xk)Trk\k-i {Xk\Zi:k-i) 


(9) 


/ gk{Zk\X)^rk\k-l{X\Z^,k-^)5X^ 

Collectively, Q and 0 are referred to as the multi-target 
Bayes filter. It was shown in |2I| that a GLMB density 


of the form 0 is closed under the Chapman-Kolmogorov 
equation (|7]) with the standard multi-target transition kernel, 
and closed under Bayes rule 0 with the standard multi¬ 
target measurement likelihood function. The GLMB is thus a 
conjugate prior for the standard multi-target tracking problem, 
facilitating the development of a closed form GLMB recursion. 

The drawback of the standard GLMB filter of | [2T| , | |2^ 
is that it does not accommodate non-standard measurement 
models, such as merged measurements, or extended targets. 
In IT), a GLMB filter for multi-target tracking in the presence 
of merged measurements was presented. Herein, we turn our 
attention to adapting the GLMB approach to tracking multiple 
extended targets. In the following section, we develop an RFS- 
based likelihood model for this problem, and then proceed to 
develop a GLMB filter for extended targets using this model. 

III. Labelled RFS -based Extended Target 
Tracking 

In this section we propose two algorithms for tracking 
multiple extended targets, based on labelled random finite sets. 
The following subsections describe the prerequisites for the 
development of the algorithms, i.e. an observation model for 
multiple extended targets, and a state space model for a single 
extended target. Based on these models, we then propose a 
GLMB filter for tracking multiple extended targets in clutter, 
as well as a cheaper approximation based on the LMB filter. 


A. Observation Model for Multiple Extended Targets 

Let us denote the labelled RFS of extended targets that exist 
at the observation time as X = {(■Ci) ^i) j ■ • •, (C|X|) ^|X|) }■ 
We formulate a measurement model based on the following 
three assumptions: 

Al. A particular extended target with state (^, 1) may 
be detected with probability pd{^,1), or misdetected with 
probability go (^, 1) = 1 - po (C, 0- 

A2. If detected, an extended target with state (^, 1) generates 
a set of detections W with likelihood g (kF|^, 1), which is 
independent of all other targets. 

A3. The sensor generates a Poisson RFS K of false observa¬ 
tions with intensity function k (•), which is independent of the 
target generated observations (i.e. K is distributed according 
to gc (K) = 

Denote hy Vi (Z) the set of all partitions that divide a finite 
measurement set Z into exactly i groups, and by U (Z) G 
Vi (Z) a particular partition of Z. For a given multi-target state 
X, denote by 0 {U (Z)) the space of association mappings 
e : C{X) ^ {0,1,..., \U (Z)\} such that 9{l) = 9{l') > 0 
implies I = V. Finally, denote by (Z) the element of the 
partition Id (Z) corresponding to label I under the mapping 9. 

Proposition 4. Under assumptions Al, A2 and A3, the mea¬ 
surement likelihood function is given by 

\x\+i 

g{Z\X)=gc{Z) Y. Y [Mzii-.O)]'' ( 10 ) 

i=l U(Z)&Vi{Z) 

e^e{u{z)) 
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where 

^U{Z) 


——’ 0(0 >o_ 

qoi^J), 0(0 = 0 


Proof: Let us first consider the case of no false detections 
(i.e. all measurements are target generated). By assumptions 
A1 and A2, the likelihood of observing a set of detections Y, 
given a set X of extended targets is given by p7) 


gD(Y\X)= Y. 9'iWi\^uh)---9'{Wix\\^ixi,lixi), 
|+)lf I Wi=Y 

( 12 ) 


where 


9'{W\^,l)x 


qoi^J), VL = 0 

PD{^,l)9{W\i,l), W^%' 


(13) 


A partition of an arbitrary set S is defined to be a disjoint 
collection of non-empty subsets of S, such that their union 
is equal to S. Note that in (12 1 , the sets kLi,..., VL|x| 
be either empty or non-empty, thus, they do not satisfy the 
definition of a partition of Y. However, the non-empty sets 
in Wi,...,VL|x| do constitute a partition of Y, hence by 
separating ( [T^ into products over the empty and non-empty 
Wj ’s, we can then write 


1^1 

gDiY\X) = [qDfY E E ( 14 ) 

i=i u{Y)ePi{Y) i<jiY-Y],<\x\ 

A PD {^jk : hk) 9 i^k {Y) ICjfc , Ij^ ) 


where Uk (Y) denotes the k-th group in partition U (Y). 
Following a similar reasoning to | p5j pp. 420], this can be 
expressed as 


gn{Y\X) = [qnfY E E 


1=1 u(Y)^Vi(Y) ee 0 (u(Y)) 


(15) 


n 

j:e(j)>o 


Pd h) 9 {^s{j) {D) \^j, 
qD {fi ) ) 


Let us now consider the case where false observations 
may also be present. By assumption A3, the set K of false 
observations has distribution gc (AT), and the sets Y and K 
are independent. The overall measurement set is Z = L U K, 
thus Z is distributed according to the convolution 


g{Z\X)= Y gc{Z-W)gD{W\X) 


wcz 


1^1 


= of ^ ^ 


wcz 


i=i u{w)ePi{w) 
ec0(u{w)) 


n 


PD (Ci4j) g i^eu) jW) 
qD (Cj I ) 


wcz i=l U{W)CVi(W) 
0C0(U(W)) 


TT Pd (^j4j) g (^e(j) (1^) 


( 16 ) 


where the last line follows from the fact that = 

nj.e(j)>o since U [W) is a partition of W. Fi¬ 

nally, this can be simplified by treating the set Z — W as 
an additional element that we append to each [W), thereby 
transforming it into a partition of Z. In doing so, the double 
summation over IL C Z and partitions U {W) G Pi (IV) up 
to size |X|, can be expressed as a summation over partitions 
of Z up to size |X| -I- 1 as follows 

l^l+i 

g(ZIX)=gc(Z)[qn]^ Y E 

i=i u(z)ePi(z) 

9ce(u(z)) 

n PD (^j4j) g (z) 
j-.eU)>o qD [ZjPjJwi 

Observe that the goi^j^j) in the denominator cancels out 
the corresponding term in the product [qo]^ when 9 (j) > 0, 
leaving one qu term for each j : 9(f) = 0. Hence, 

( [TT] ) can be equivalently expressed in the form ([T 0 li-(rT]i. ■ 

In general, an exact calculation of the likelihood ([g will 
be numerically intractable, because the sets of measurement 
partitions and group-to-target mappings can become extremely 
large. However, it has been shown that in many practical 
situations, it is only necessary to consider a small subset 
of these partitions to achieve good performance |[T 6 ). 

Additionally, the set of group-to-target mappings can be 
substantially reduced using a ranked assignment algorithm, 
thereby cutting down the number of insignificant terms in the 
likelihood even further. 


B. Extended Target State-space Model 


In this section we describe the extended target state space, 
and the class of probability distributions used to model a single 
extended target. We begin by introducing some notation; 

• IR+ is the space of positive real numbers 

• K" is the space of real n-dimensional vectors 

• is the space of n x n positive definite matrices 

• S" is the space of n x n positive semi-definite matrices 

• QAM. ( 7 ; a, fS) is the gamma probability density function 
(pdf) defined on 7 > 0 , with shape a > 0 , and inverse 
scale /3 > 0; 

gAM(T,aM) = 

T(a) 

• Af (x;m, P) is the multivariate Gaussian pdf defined on 
X £ K", with mean m £ K" and covariance P £ S" 


Af (x; m, P) = 


V(‘2^r \p\ 
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y) is the inverse Wishart distribution defined 
on X S with degrees of freedom v > 2d, and scale 
matrix V S 


iWd {x;v,v) = 


i^r 






where (•) is the multivariate gamma function, and tr (•) 
takes the trace of a matrix. 

• Id is the idenity matrix of dimension d. 

• A(^ B is the Kronecker product of matrices A and B 
The goal is to estimate three pieces of information about each 
target; the average number of measurements it generates, the 
kinematic state, and the extent. We thus model the extended 
target state as the triple 


e=(7,a;,x)GK+xM"^ xSt+, (18) 


where 7 G K+ is the rate parameter of a Poisson distribution 
that models the number of measurements generated by the 
target, x G M"* is a vector that describes the state of the 
target centroid, and x G S++ is a covariance matrix that 
describes the target extent around the centroid. The density 
of the rate parameter is modelled as a Gamma distribution, 
the kinematics as a Gaussian distribution, and the covariance 
of the extent as an inverse-Wishart distribution. The density 
of the extended target state is thus the product of these three 
distributions, denoted as a gamma Gaussian inverse Wishart 
(GGIW) distribution on the space K+ x x given by 


P{0 = P il) P {x\x) P ix) 

= GAM ( 7 ; a,l3) X M {x] m,P ®x) x IWd (x; v, V) 
= GGTW{i-X) (19) 


where C, = {a, /3, m, P, v, V) is an array that encapsulates the 
GGIW density parameters. We now describe the prediction and 
Bayes update procedures for a GGIW distribution representing 
a single extended target. 

1) Prediction : The predicted density p_|_ (•) of an extended 
target is given by the following Champan-Kolmogorov equa¬ 
tion 

p+{o = I fimpiOde, (20) 

where p{-) = GGHW {-X') 1^ th® posterior density at the 
current time with parameters = (a', /3', m', P' , v' , V), and 
/(•j) is the transition density from the current time to the 
next time. This has no closed form solution, hence we resort 
to making a GGIW approximation for (^). We start by 
assuming that the transition density can be written as the 
product E) 

/ I?') = /t ( 7 I 7 ') fx {x\x, x') fx (xlxO , (21) 

which yields the following predicted density 

P+ iO = j GAM ( 7 ; a', X) fj ( 7 I 7 ') dX 

X J Mix';m',P'i^x) fx{x\x,x')dx' 

X J IWd (x'; v', V) fx (xlx') dX- (22) 


If the dynamic model is linear Gaussian with the form 
fx {x\x, x') = N (x; {F ® Id) x', Q (S> x)’ the kinematic com¬ 
ponent (i.e. the second line in ( [ 22 ] i) can be solved in closed 
form as follows 

J Af{x'; to', P' (g) x)fx {x\x, x') dx' = M {x;m, P x), 
m={F® Id) to', P = FP'F^ + Q. (23) 


However, closed forms still cannot be obtained for the mea¬ 
surement rate and target extension components, which can 
be addressed by the use of some additional approximations. 
For the measurement rate component we use the following 
approximation proposed in 113^ 

J GAM ( 7 '; a', X) /t ( 7 I 7 ') dxk-i « GAM ( 7 ; a, j3 ), 


a = 


P 



(24) 


In the above, p = is an exponential forgetting factor 

with window length w > 1. This approximation is based on 
the heuristic assumption that £’[ 7 ] = £[ 7 ^]^ Var( 7 ) = 
Var( 7 ') X p, i.e. the prediction operation retains the expected 
value of the density, and increases its variance by a factor of 
P- 

For the extension component we use the following approx¬ 
imation, as proposed in Q, 

J IWd (x; v', V) f^ (xlxO dx' - IWd (x; , 


V = V = ] v'. (25) 

v—d—1 

Similary to the measurement rate, this approximation assumes 
that the prediction retains the expected value and reduces the 
precision of the density. For an inverse-Wishart distribution, 
the degrees of freedom parameter is related to the precision, 
with lower values yielding less precise densities. A temporal 
decay constant r is thus used in ( |25| ) to govern the reduction 
in the degrees of freedom. Based on the calculated value for v, 
the expression for V retains the expected value of the inverse- 
Wishart distribution through the prediction. 

The above yields an approximate representation of the 
predicted GGIW density X) « GGIW {f; f), where f = 
(a, P, TO, P, V, V) is the array of predicted parameters defined 
by equations and 

2) Update: In the proposed GGIW-(G)LMB filter, each 
extended target will need to undergo measurement updates 
using various subsets of the measurement received on each 
scan. In what follows, we describe the update procedure for a 
single target with predicted density p(-) = GGIW {-X)^ for 
a given extended target generated measurement set W. The 
first step is to calculate the mean and scale matrix of W, 
the innovation, innovation factor, innovation matrix and gain 
vector; 

' ' wew 

T* = {w — w) {w — w)'^ , 

wGW 


(27) 
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e = w — {H ® Id) m, (28) 

* = + WY 

N = (30) 

K = PH'^S-\ (31) 


The posterior GGIW parameters are then given by C,w = 


{aw, /3w,'mw, Pw,vw,Vw), where 

aw = a+\W\, (32) 

/?ty=/3 + l, (33) 

mw = m + {K ® Id) e, (34) 

Pw = P - KSK^, (35) 

v^=v + \W\, (36) 

Vw = V + N + ^!. (37) 


The GGIW-(G)LMB filter also requires calculation of the 
Bayes evidence for each single-target update, as they are 
needed when computing the weights of the posterior GLMB 
components. This is given by the product of the following two 
terms. 


77 {W',(,Cw) 

Vx,x {^, C, Cw ) 


1 r(ai^)/3“ 
\W\lT{a)(3^-'^ 

si \Vwf^Td{^) 


(38) 

(39) 


Note that the measurement rate component ( [38] l corresponds 
to a negative-binomial pdf, and the kinematics-extension com¬ 
ponent ( |39l l is proportional to a matrix variate generalized beta 
type n pdf | [T3) . 


C. GLMB Filter for Extended Targets 


We now present a GLMB filter for extended targets, based 
on the measurement likelihood and state space models de¬ 
scribed in the previous sections. The GLMB filter consists 
of two steps, prediction and update. Since we are using the 
standard birth/death model for the multi-target dynamics, the 
prediction step is identical to that of the standard GLMB 
filter derived in ID- For completeness, we shall revisit the 
final prediction equations, and the reader is referred to pT| 
for more details. Denote by pg (^, 1) the probability that a 
target with state {^,1) survives to the next time step, and by 
9s (C) 1) = ^ — PS (^) 0 '^he probability that a target does not 
survive. The birth density is an LMB with label space B, 
weight wb {■) and single target densities pB (-, 0- '^he multi¬ 
target posterior is a GLMB of the form Q with label space 
L, then the predicted multi-target density at the next time step 
is the GLMB with label space L+ = L U B given by 


7r+ (X) = A(X)^u;« (/:(X)) 


cGC 


P+’ (•) 


n X 


(40) 


where 


(L) = Wb (L — L) (L n L), (41) 

P^:^ (?, l) = 1l (0 pf (e, 0 + (1 - 1 l {l))PB (?, l), (42) 


Ps^ (C ,0 = 


!ps {i',i)di' 


pf (() 


pf{i)= I 


' ;i IP 


(J)= 


(c) 


n J 


Ps 




/CL 


4 "^ ( 0 = J Qs {^,l)p^^^ ii,l)dC 


(43) 

(44) 

(45) 

(46) 


The function f{-\-,l) is the single-target transition kernel, 
which in this case is the GGIW transition defined in Section 
IIII-BII 

Clearly, the difference between the standard GLMB and 
extended target GLMB filters will lie in the measurement 
update procedure, since the measurement likelihood function 
has a different form. The update for the extended target GLMB 
is given by the following proposition. 

Proposition 5. If the prior is a GLMB of the form (1^ then 
under the extended multi-target likelihood function ( (7Q| ), the 
posterior is a GLMB with label space L+ = L U B, given by 


\x\+i 

TT {x\z) = a{x)J2J2 E ^41) (^)) 

cGC i=l U(Z)eVi{Z) 

9eeiu(z)) 


X 




(47) 


where 

{c,9) 

^U(Z) 


{L) = 


(L) 


(c,0) 

Pu{z) 


n L 


LI + 1 r r 

LEE E {J) 


cGCJCL i=l U{Z)€:Vi{Z) 


Pu(z) 


{U\U{Z)) = 


MS) 


9g0(U(Z)) 

(48) 

{U)Mz) {U-.d) 

dufz) (0 

(49) 

{^,i)Mz) ii,i-,d)d^ 

(50) 


in which ipuiz) (Ci d) is given by {11 L 

Proof: The product of the prior distribution and the 
likelihood is 


rT{X)g{Z\X) 


\x\+i 

= A(X)pc(^)EE E 

cGC i=l U{Z)eVi{Z) 
9&e{u(z)) 
1 X 


X 


P^^"’ (OV'WCX) {■■,d) 


l^l+i 


= A(x)gc(^)EE E 


cGC i=l U[Z)^Vi(Z) 
9&0iU(Z)) 

P^^’<^H-\U{Z))pI%] {■)]'' 
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|£(X)| + 1 

= A(x)gc(^)E E E 

cGC i=l U{Z)eVi(Z) 

ee 0 {u(z)) 


^ufz) (•) 


(Z)) . (51) 


In what follows, we use the simplifying notation 

(^Xj , )), ^i:j — (^15 ■ ■ • 7 ) and x\-j — 

{xi, ... ,Xj) to denote vectors, with |(a;, and {/i:j} 

denoting the corresponding sets. The set integral of ( [5T| l with 
respect to X is given by 

17T{X)g{Z\X)6X 

l.\C{X)\+l 

= gc{Z)Y, E E ^iX)w<-^Hc{X)) 


1) LMB Prediction: If the multi-target density at the 
current time is an LMB of the form 0 with parameters 
I 1^^^, and the multi-target birth model is an LMB 

with parameters | | ’ then the predicted multi¬ 

target density at the next time step is an LMB with parameters 

I () I , with L+ = L U B comprising both surviv- 
V \ / J 

ing and birth components 

= {(’■A'Pjy },„ u{(r«,p«)} (53) 


where 

r+]s = VS 

(0 IPs{C,l)fiiW,l)/HC)dC 


P+\s = 


VS (0 


(54) 

(55) 


cGC ‘ 


i=i u{z)eVi{z) 
0G0(UiZ)) 


vs (1)= I I PS (C', 1) f (eir, 1) P^'^ iC, 1) d^d^- (56) 


G,s)' 

C{X) 

Vu{z) 



n X 


sx 


= 9c 


p^^’0H-\U (Z)) 
i +1 

(z)EEi EE E 


ceCj=0 ■' i=lU{Z)&Vi(Z) 

e&e(u(z)) 




(c,0)' 


Vu{z) 

Ji 


i-lUiZ)) 

{c 




d (a;i 


= 5c(^)EE E E 

cGCLCL i=l U{Z)^Vi{Z) 
0 G 0 (U{Z)) 


(cfi) 

Vu(Z) 


1 L 


■ (52) 


In the above, the second line is obtained by applying 
proposition 2 from pT[ , and taking the parts that are only 
label-dependent outside the resulting integral. The last line is 
obtained by observing the fact that the distinct label indicator 
function limits the summation over j : 0 —>■ c» and li-j G L-^ 


to a summation over the subsets of L. Substituting (51 1 and 
into yields the posterior density (gT}. ■ 


Note that this result establishes that the GLMB is a conju¬ 
gate prior with respect to the extended multi-target measure¬ 
ment likelihood function. 


That is, to obtain the predicted LMB, we simply take the 
union of the predicted surviving tracks and the birth tracks. 
This is much cheaper to compute than the GLMB prediction, 
since it does not involve the sum over subsets of L which 
appears in ( |45| ). The reader is referred to for more details. 

2) LMB to GLMB Conversion: The update step requires 

. ^converting the LMB representation of the predicted multi¬ 
target density to a GLMB representation. The predicted 
LMB 7r+ = can be converted to a single 

component GLMB, given by In principle, this involves 
calculating the GLMB weight for all subsets of L+, however, 
in practice, approximations can be used to reduce the number 
of components and improve the efficiency of the conversion. 
These include methods such as target grouping, truncation 
with k-shortest paths, or sampling procedures. For more details 
see p^ , which expresses the converted density in (5-GLMB 
form, essentially enumerating all possible subsets of the tracks 
appearing in the predicted LMB. 

3) Approximating GLMB as LMB: After the measure¬ 
ment update step, the posterior GLMB (which is given by 
Proposition can approximated by an LMB with matching 
probability hypothesis density, with parameters 


D. LMB Filter for Extended Targets 

The key principle of the labelled multi-Bernoulli (LMB) fil¬ 
ter is to simplify the representation of the multi-target density 
after each update cycle, in order to reduce the algorithm’s com¬ 
putational complexity. Instead of maintaining the full GLMB 
representation from one iteration to the next, we approximate 
it as an LMB representation after each measurement update 
step. In the subsequent iteration, we carry out the prediction 
step using this LMB representation, before converting the 
predicted LMB back to a GLMB in preparation for the 
next measurement update. Thus there are three modifications 
needed to turn the GLMB filter into an LMB filter; 1) replace 
the GLMB prediction with an LMB prediction, 2) convert the 
LMB into a GLMB representation in preparation of the update, 
and 3) approximate the updated GLMB density in LMB form. 


"^'^=EE E (57) 

cGC i=l U(Z)^Vi(Z) 

e&0{u{z)} 

, \L\+l 

(^) = ^(0 E E E ’^uS) (^) (0 (c, 0 • 

cGC i=l U(Z)^Vi{Z) 

e&e(u{z)) 

(58) 

The existence probability corresponding to each label is the 
sum of the weights of the GLMB components that include 
that label, and its pdf becomes the weighted sum of the 
corresponding pdfs from the GLMB. Thus, the pdf of each 
track in the LMB becomes a mixture of GGIW densities, 
where each mixture component corresponds to a different 
measurement association history. To avoid the number of 
components growing too large, it is necessary to reduce this 


















mixture by a process of pruning and merging. This can 
be carried out using the techniques proposed in Og, |T7), 
which have been previously applied in the context of mixture 
reduction for the GGIW-CPHD filter |T4}. 


IV. Implementation 


This section provides more details (including pseudo-code) 
on the implementation of both the GLMB and LMB extended 
target filters. We begin in Section |IV-A| by describing the 
prediction and update steps of the GGIW-GLMB filter, then 
in Section |IV-B we describe the modifications necessary to 
implement the GGIW-LMB filter. To simplify the presenta¬ 
tion of the pseudo-code, the following functions are used to 
encapsulate some of the lower-level procedures; 

• Poisson {y, A): Poisson pdf with mean A computed at 
each element of the array y. 

• Allocate (n, w): Randomised proportional allocation of 
a scalar number of objects n, into a fixed number of bins 
with weights given by the array w. 

• Normalise (w): From a set of unnormalised posterior 
component weights w, compute the posterior cardinality 
distribution and normalised components weights. 

• PredictGGlW (Af); Predict the GGIW X up to the cur¬ 
rent time, using the method in Section |III-B1| 

• UpdateGGlW (Y, IF): Update the prior GGIW Y with 
measurement set W, using ([26]i-(37 1 . 


Note that throughout the pseudo-code, $ and fl are used to 
denote GLMB densities, and $ and fl denote LMB densities. 
A GLMB is represented as a data structure containing four 
arrays; $.Ar contains the single target pdfs, ^.L contains the 
target labels, <i).w contains the component weights, and ^.p 
is the cardinality distribution. Up to three indices are used 
to identify elements within these arrays; the first indicates 
cardinality, second is the component index, and third is the 
target index. For example, and are the 

pdf and label of the f-th target in the m-th component of 
cardinality n, is the weight of the m-th component 

of cardinality n, and is the value of the cardinality 

distribution corresponding to n targets. 

An LMB is represented as a structure containing the arrays 
$.Ar (single target pdfs), $.L (target labels), and $.r (target 
existence probabilities). Since the LMB does not have multiple 
components, a single index is sufficient to identify the pdf, 
label and existence probability of any particular target. 


A. GGIW-GLMB Filter 

1) Prediction: To compute the predicted GLMB density, 
we predict the individual target pdfs forward using (|22)l-([25]l. 
For each component c in the previous density, the values from 
(|44)l and (46 1 are used to construct the following cost matrix 


/ 


= - log 


ih) 


V i^n) Qs’ (In) y 


(C 

Ys 


ih) \ 


(59) 


which is denoted in the pseudo-code by 
CostMatrixSurvive (Af,ps). The cost matrix is used 


to calculate components for the predicted GLMB of surviving 
targets. This is done by constructing a directed graph, where 
every element of the cost matrix becomes a node. Each 
node has two outgoing edges, one ending at each node 
in the following row. We then generate the K shortest 
paths from the top row to the bottom row, denoted by the 
function ShortestPaths (C, Ff). Each path corresponds to 
a predicted component, that comprises the targets associated 
with the rows in which the first column was visited. 

After generating the surviving target density, it is multipled 
by the GLMB of spontaneous births, yielding the overall 
prediction. Pseudo-code for the GLMB prediction is given in 
Eigure 


U = PredictGLMB ($, <1>b, N,ps) 


2 

3 

4 

5 

6 

7 

8 
g 

10 

11 

12 

13 

14 

15 


A = MecUl (^.p); 

— Allocate {N, Poisson (1 : ; A)); 

k = Ones (1, 
for n — 1 : do 

'^n,l:|<I..x(»)| = Allocate 
for i = 1 : do 

C= CostMatrlxSurvive ; 

P = ShortestPaths (C, Kn^i); 
for (p, c) E P do 

« = U; 1 < j < n,p{j) = 1}; 
m = |s|; 

^ PredlctGGIW ) ; 

k (m) = k [m) + 1; 


16 ^s-P) — Normalise (w); 

17 Q = $5 X 


Figure 1. Pseudo-code for GGIW-GLMB prediction. Inputs: $ is the 
posterior GLMB at the previous time step, is the GLMB density of 
spontaneous birth targets, N is the number of prediction components to 
generate, and ps is the target survival probability. Output: is the predicted 
GLMB density at the current time. 


2) Update: Similarly to the extended target (C)PHD filter, 
the main barrier to implementing the extended target GLMB 
filter is the fact that the posterior density in (47i involves a 
sum over all partitions of the measurement set. Even for small 
measurement sets, exhaustively enumerating the partitions is 
usually intractable, because the number of possibilities (given 
by the Bell number) grows combinatorially with the number 
of elements. Therefore, to make the filter computationally 
tractable, the first step in the update procedure is to reduce the 
number of partitions to a more managable level, by removing 
those that are infeasible. 

Ideally, the retained partitions should be those that give 
rise to GLMB components with the highest posterior weights, 
such that the effect of truncation error is minimised. Although 
it is difficult to establish a method that can guarantee this, 
the use of clustering techniques to generate the most likely 
partitions has been shown to produce favourable results |T3), 

1161. In our implementation of the GGIW-GLMB filter, we use 
a combination of distance-based clustering and the expecation- 
maximisation algorithm to generate a set of feasible partitions 
of the measurements, in a similar manner to m and Hi- 
Eor the pseudo-code in Eigure this is encapsulated by the 
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function FeasiblePartitions (Z). 

After generating the feasible partitions, each unique group¬ 
ing of measurements is used to update the single-target GGIW 
pdfs in the GLMB density, using The implementa¬ 

tion is simplihed by assuming that the detection probability 
is dependent on the target label only (i.e. po (^, 1) = pu (1)), 
yielding the following closed form expression for ([50|l, 


(0 = 


Pd (0 {Z) (Z) 




1){Z) 


(60) 

where rjj (•) and rjx^x (') given by ( [38| l and (39i. 

For each predicted component c, a cost matrix is calculated 
for the assignment of measurement groups to targets. This 
matrix is of the form M(°(] , where is a 

|X| X \IA {Z)\ matrix and is a |X| x |X| matrix, the 

elements of which are 


= 
* J 


■ log (vu%) ih) 


(c) _ j - log {qd ik)), i= j 

^ j 


^^3 = 


( 61 ) 

(62) 


where denotes the Bayes evidence (601, from the 


' ' ' I_I 

update of the target labelled U within component c, using 
the measurement group Uj (Z). Note that each row in 
corresponds to a target, and each column corresponds to 
either a group of measurements in Z or a misdetection. In 
the pseudo-code, the calculation is denoted by the 

function CostMatrixAssign (if (Z) ,X,pd)- 

Murty’s algorithm is then used to generate highly weighted 
assignments of measurement groups to targets, denoted by 
the function Murty (C, n), which generates the n-best ranked 
assignments based on the cost matrix C. Each assignment 
returned by Murty’s algorithm forms a component in the pos¬ 
terior GLMB density. Pseudo-code for the update procedure 
is given in Figure 

3) Track Extraction and Pruning: After the update, labelled 
target estimates are extracted from the posterior GLMB, which 
are used to update a table of reported tracks. This is done by 
hnding the maximum a-posteriori estimate of the cardinality, 
then selecting the highest weighted GLMB component with 
that cardinality. For those labels in the selected GLMB com¬ 
ponent that are already present in the reported track table, the 
current estimates are appended to the corresponding reported 
tracks. For labels that are not present, new tracks with those 
labels are inserted into the reported track table. The posterior 
GLMB is then pruned by retaining the top M components 
with highest weights. More details and pseudo-code for these 
procedures can be found in |[T]. 


B. GGIW-LMB Filter 

1) Adaptive Birth: In previous work p4) we used a static 
model for target birth, which meant that new targets could 
only be initiated around pre-determined locations. To alleviate 
this restriction, we now use an adaptive target birth model that 
allows for new targets to appear anywhere in the state space. 
The adaptive birth density for the GGIW-LMB hlter closely 


n = UpdateGLMB (3', Z, N,pd) 


2 

3 

4 

5 


r 


10 


A = Mean 

— Allocate {A^, Poisson (1 : |$.A'| ; A)); 

V = FeasiblePartitions(Z); 
for n — 1 : Ifi.Xj do 
k^l-, 

= Allocate 

for 2 — 1 : I I do 

^i:\V\ — Allocate (Ji,Unif(l : |P|)); 
for — 1 : |P| do 

C — CostMatrixAssign (Z) , ; 


11 

12 

13 

14 

15 

16 


A = Murty {C, Kj); 
for (a, c) 6 X do 

for I = 1 : n do 

if a{l) > 0 then 

Qj^(n,fc,0 


17 


= UpdateGGIW 






U) 

a{l) 



18 

19 

20 
21 


else 

k — k 1; 


22 (Q.m, Q.p) = Normalise (it?); 


Figure 2. Pseudo-code for GGIW-GLMB update. Inputs: is the predicted 
GLMB density at the observation time, Z is the current measurement set, and 
N is the maximum number of posterior components to generate, pp) is the 
detection probability. Output: O is the posterior GLMB density. 


resembles the adaptive birth density of the standard LMB hlter 
i.e. measurement clusters that are far away from any 
existing tracks are likely to correspond to new born targets, 
while clusters in the proximity of existing tracks are likely to 
have originated from these tracks. Hence, the adaptive birth 
density is dominated by the clusters corresponding to possible 
new born targets. The reader is referred to p3| for additional 
details about the adaptive LMB birth density. 

2) Prediction: The LMB prediction involves computing the 
predicted GGIW for each target in the density, and multiplying 
the existence probability of each target by its probability 
of survival. This is signihcantly cheaper than the GLMB 
prediction, since it does not require generation of components 
using k-shortest paths. The predicted LMB density is then 
given by the union between the surviving and birth densities, 
where the birth density may be static or adaptive depending 
on the scenario requirements. The pseudo-code for the LMB 
prediction is given in Figure 


n = PredictLMB (#, isiPs) 

1 for i = 1 : do 

2 I's.L'*' = 

3 = PredictGGIW ($.X(’)) ; 

4 I's.rt’) =p, (•i.L(’)) X i.rt’); 

5 0 = 


Figure 3. Pseudo-code of GGIW-LMB prediction. Inputs: $ is the posterior 
LMB at the previous time step, is the LMB of spontaneous births, and 
PS is the target survival probability. Output: O is the predicted LMB at the 
current time. 
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3) Update: The LMB update consists of splitting up the 
predicted LMB into groups of well separated targets. This is 
done based on the current measurement set, and the obser¬ 
vation that a pair of tracks can be considered as members 
of different independent clusters if no single measurement 
falls inside the gating region of both tracks simultaneously. 
The clustering procedure is denoted in the pseudo-code as 
ClusterTracks ($, Z), and more details on how this is 
performed can be found in 

Following clustering, the predicted LMB for each group is 
converted to a GLMB using (|^-Q. For an LMB density l> 
and a given group of labels L, this conversion is denoted by the 
function LMBtoGLMB Each of these is updated using 

the UpdateGLMB procedure in Figure]^ The posterior GLMB 
for each group of targets is then approximated in the form of an 
LMB using (|^-@ followed by a mixture reduction step for 
each track (denoted by the function ApproximateLMB ($)). 
Finally, the union of the posterior LMBs is taken across all 
target groups to obtain the overall posterior LMB density. The 
pseudo-code for the LMB update is given in Figure 


n = UpdateLMB Z, N,pd^ 

■L C = ClusterTracks (<i*, ; 

2 0 = 0 ; 

3 for G C do 

4 = LMBtoGLMB (#, W) ; 

5 Ovv^ = UpdateGLMB (^>Wi 

6 Ow^ = ApproximateLMB (Opv^); 

7 O = O U 0^/^/^ 


Figure 4. Pseudo-code for GGIW-LMB update. Inputs: is the predicted 
LMB density at the current time, Z is the current measurement set, and N is 
the maximum number of posterior components to generate in the GLMB 
update, is the detection probability. Output: O is the posterior LMB 
density. 


Note that the computational saving of the LMB filter 
depends on the assumption that the number of targets (and 
hence the number of GLMB components) in each group will 
be relatively small. In this case, the total number of GLMB 
components that need to be processed across all target groups 
will be significantly lower compared to the full GLMB filter. 


V. Simulation Results 


In this section, the performance of the GGIW-GLMB, 
GGIW-LMB, and GGIW-LMB with adaptive birth process 
(GGIW-LMB-ab) is compared to an extended target CPHD 
filter | |T4| using the cardinality estimation error and the opti¬ 
mal sub-pattern assignment (OSPA) distance |30|. Since the 
standard OSPA only penalizes cardinality and state errors, 
a modified version of the OSPA metric m Section VI] 
incorporating measurement rates and target extent is used in 
the evaluation. 

The targets follow the dynamic model 


Xk+i = {Fk+i\k ® Id) Xk + Vk+i (63) 

where Vk+i Af (O, Qfc-i-i|fc) is a d x 1 independent and 
identically distributed (i.i.d.) process noise vector. Id is the 


identity matrix of dimension d, and Fk+i\k and Qk+i\k are 



■ 1 T 

1 'Y'2 

2^ 



0 0 

T 

7 


. 0 0 

e-T/e 


Qk-\-l\k 

E^ (l- 

e-2T/e^ 

diag([ 0 0 1 ]) <ZXk+i 


In the above, T is the sampling period, E is the scalar accel¬ 
eration standard deviation, and 6 is the manoeuvre correlation 
time. In these simulations, we use parameter values of T = 1 s, 
0 = 1 s and S = 0.1 raja^. 

The forgetting factor used by the filters in (24i for the 
prediction of target measurement rates is set to p, = 1.25, and 
the temporal decay constant in ( |25| ) for the prediction of the 
target extension is r = 5 s. The probability of target survival 
is set to ps = 0.99. For the GLMB and LMB with static 
birth, the parameters of the gamma components are ao = 10 
and Po — 1, the inverse-Wishart component parameters are 
j/Q = 10 and Vo = 100 x G, and the kinematic components 
have means mo which are located close to the true target 
starting positions with zero initial velocity and covariance 
Pq — diag([ 10 2.5 ]) . The same values for ao, /3o, i^o 

and Pq are used in the LMB filter with adaptive birth, however, 
the values of mg and Vq are computed on-line. 

The measurement model for a single detection is 


Zk = {H ® Id) Xk + Wk (64) 

where H = [ 1 0 0 ], and Wk ~ Af {0,Xk) is i.i.d. 

Gaussian measurement noise with covariance given by the 
target extent matrix Xk- If detected, a target generates a 
number of measurements from the model ( |64l i, where the 
number follows a Poisson distribution, the mean of which 
may be set differently for each target. In addition, clutter 
measurements are simulated as being uniformly distributed 
across the surveillance region, where the number of clutter 
points is Poisson distributed with a fixed mean. 

Three scenarios were simulated, the first two of which were 
used in | |T4) to compare the performance of the GGIW-PHD 
filter and the GGIW-CPHD filter (note that all scenarios are 
2-dimensional, i.e. d = 2). Scenario 1 runs for 200 time steps, 
and consists of four targets that appear/disappear at different 
times. The targets generate measurements with a detection 
probability of pu = 0.8 and the clutter measurements follow 
a Poisson distribution with a mean number of 30 per time 
step. Due to lower detection probability, higher clutter rate, 
and target birth/death, the estimation of the cardinality is 
challenging in this scenario. Scenario 2 runs for 100 time 
steps and consists of two targets that are present for the 
entire scenario. The two targets are spatially well separated 
at the beginning, then move in parallel at close distance, 
before separating again towards the end. In this scenario, the 
detection probability is pu = 0.98 and a mean value of 10 
Poisson distributed clutter measurements occur at each time 
step. This scenario is used to illustrate the filter performance 
for the difficult problem of tracking closely spaced targets 
|T3|, fig. Since the target-generated measurements are close 
together, they often appear as a single cluster in the sensor 
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data, rather than multiple separate clusters. Figure depicts 
the true trajectories of the targets for both scenarios. 



Figure 5. Simulated true target tracks. In scenario 1 (left) all tracks start in 
the origin. In scenario 2 (right) the tracks stait on the left. 

Scenario 3 is used to test the so-called ‘spooky’ effect m- 
The scenario has two targets that are spatially separated by at 
least 1km for all 50 time steps. The probability of detection 
was set to pjj = 0.9 and clutter Poisson rate was 10. The 
measurements were generated such that one target is always 
detected, and the other target is detected on all time steps, 
except for steps 20, 40 and 41. 

For the first two scenarios, 1000 Monte Carlo runs were 
carried out in order to compare the performance of the four 
different GGIW filters: GLMB, LMB, LMB with adaptive 
birth process (LMB-ab), and CPHD. The GGIW-PHD filter 
is omitted from the comparison because previous work has 
shown that both the CPHD and GLMB filter outperform the 
PHD filter p4| . Figure shows the mean OSPA distances, 
and the mean cardinality errors for scenario 1. The GLMB 
and LMB filters have approximately equal performance. The 
LMB-ab filter has slower convergence due to the unknown 
birth density, however the filter eventually reaches the same 
error as the GLMB and LMB with known birth density. This 
is expected since the other three filters have the advantage 
of knowing the region where new targets will appear. The 
CPHD filter can match the GLMB and LMB in terms of the 
cardinality error, however the mean OSPA is larger. 

The execution times for our Matlab implementation of the 
algorithms (mean + one standard deviation) for Scenario 1 are 
3.95 + 3.41s for the GLMB filter, 0.19 + 0.29s for the LMB 
filter, and 2.20 + 0.47s for the CPHD filter. Since the LMB-ab 
filter only uses clusters with more than four measurements as 
birth candidates, it is even faster than the LMB filter. Since 
the LMB filter partitions the tracks and measurements into 
approximately statistically independent groups ||3T), its 
computation times are less than those of the CPHD filter. 

Figure [7] shows the mean OSPA distances and mean cardi¬ 
nality errors for scenario 2. Similar to scenario 1, the GLMB 
filter slightly outperforms the LMB and CPHD filter. Again, 
the cardinality estimate of the LMB-ab filter takes longer to 
converge to the correct value, however, the LMB-ab filter has 
lower OSPA distance and smaller cardinality error than the 
LMB filter after time 78. This is due to the fact that, in some 
of the runs, the LMB filter lost one of the tracks because 
the measurement clusters were very close together. Even after 
the targets move apart, the filter with static birth density is 
unable to start a new track on the lost target. However, the 
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Figure 6. Cardinality error and OSPA metric for scenario 1 (mean values of 
the 1000 Monte Cai'lo runs) 


filter with adaptive birth density has the capability to start 
a new track at the lost target’s current location, leading to 
improved performance. 



Time 



Time 


Figure 7. Cardinality error and OSPA metric for scenario 2 (mean values of 
the 1000 Monte Carlo runs) 


In scenario 3 we compare the PHD, CPHD and LMB 
filters. The estimated weights (for the PHD and CPHD) and 
existence probabilities (for the LMB) for a single run are 
shown in Figure The PHD filter suffers from a positive 
bias (weight around 1.1), and the weight drops quickly when 
there are missed detections. The CPHD clearly suffers from 
the ‘spooky’ effect p9[ , as the weight of the detected target 
increases when the other target is misdetected. In comparison, 
the LMB filter performs better, as the probability of existence 
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of the detected target is unaffected when the other target is 
not detected. Also, the decrease in the existence probability 
following the missed detections is more conservative compared 
to the decrease in the weights for the PHD and CPHD biters. 



Time 



Figure 8. Result from spooky effect scenario. The lines show the estimated 
weights for two spatially separated targets. The distance between the two 
targets is 1km, at time steps 20, 40 and 41 there are missed detections for 
target 1. 


VI. Experimental Results 

To demonstrate the proposed method on a real world sce¬ 
nario, we applied the GGIW-LMB biter to pedestrian tracking 
with laser rangebnders. In contrast to targets with distinct 
shapes, such as vehicles or buildings, pedestrians do not 
exhibit specibc structure in laser scans, appearing instead as 
a random cluster of points. Hence, the GGIW measurement 
model is well suited to this application, since it assumes 
that the measurements are distributed normally around the 
target centroid. For this experiment, two pedestrians were 
recorded while walking on a parking lot using three Ibeo 
Lux laser sensors, which are mounted in the front bumper 
of the vehicle. Before being passed to the tracking biter, the 
laser returns from each sensor were thinned by removing 
measurements that lie outside the region of interest, thus 
excluding measurements from parked vehicles. The sensor was 
stationary during the experiment, and both pedestrians were 
wandering around the surveillance region, and are in close 
proximity to each other mid-way through the scenario. Figure 
[T0| shows camera footage from this instant. 



Figure 10. Pedestrian tracking scenario: The two pedestrians are getting 
close at around t = 6.5 seconds. 


A constant velocity dynamic model is used to track the 


pedesttians, i.e., the parameters of (63 i are d = 2 and 
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where the sampling period is T = 0.08 s and the standard 
deviation of the process noise is cr = 4m/s^. Due to the 
decreased dimension of the motion model, the measurement 
matrix in (64 1 is Hk = [ 1 0 ], and since the birth locations 


are unknown, the LMB biter with adaptive birth model is used 
in this scenario. Other biter parameters are similar to those 
used in Section ( [V| ). 

The results from the GGIW-LMB-ab biter are depicted 
in Figure The dashed lines show approximate ground 
truth trajectories, which were obtained by manually labelling 
the pedestrians in the raw laser scans, and the solid lines 
show the estimated ttajectories. In addition, two-sigma ellipses 
representing the target extents as well as the corresponding 
measurements are also shown for selected time instants. 

The biter is able to track both pedestrians continuously, 
even when they are very close. Especially in this situation, 
the multi-target representation as (G)LMB facilitates bnding 
consistent association hypotheses and maintaining tracks over 
time. Note that the strongly buctuating measurements on 
different moving parts of the human body, such as legs, 
arms and torso, make precise estimation of the target centroid 
positions difficult, in both the manual labelling process, and for 
the tracker itself. This explains most of the deviation between 
labelled ground truth and estimated trajectories. The estimates 
of the pedestrian extent vary over time, as demonstrated 
by the changing ellipses in Figure This is again due to 
buctuating measurements, which can be attributed mostly to 
leg movement. When the targets are close to the sensors, and 
the laser rangebnders provide detailed scans of the legs, a 
periodic adaption of the target extent following the motion of 
the legs with each stride could be observed. 


VH. Conclusion 

In this paper we have proposed two algorithms for tracking 
multiple extended targets in clutter, namely the GGIW-GLMB 
and GGIW-LMB biter. Both are based on modelling the prob¬ 
lem using labelled random bnite sets, and gamma Gaussian 
inverse Wishart mixtures. The proposed algorithms estimate 
the number of targets, and their kinematics, extensions and 
measurement rates. The major advantage of these methods 
over the previously developed GGIW-(C)PHD biters is the 
inclusion of target labels, allowing for continuous target backs, 
which is not directly supported by the (C)PHD biters. 

Of the two proposed algorithms, the GGIW-GLMB biter 
is more accurate as it involves fewer approximations, but 
it is also more computationally demanding than the GGIW- 
LMB biter. Simulation results demonstrate that the algorithms 
outperform the GGrW-(C)PHD biters, especially in cases 
where the performance of the CPHD biter is degraded due 
to the spooky effect. Finally, we have also demonsbated that 
the biter performs well in a real-world application, in which 
laser rangebnders are used to back pedesbians. 
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Figure 9. Ground tmth (dashed) and estimated (solid) pedestrian trajectories. Estimated two-sigma ellipses and corresponding laser measurements (black) 
are plotted for selected time steps. 
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